function [t, A_lo, A_corr, c_lo, c_corr, c_num] = verify(k0, Nmodes, asy, figs)


if (nargin < 3)
    asy = 0;
end

if nargin < 2
    Nmodes = 2;
end

if nargin < 4
    close all
    figs = 1;
end


% load parameters
p = params;

t1 = logspace(-6, 0, 800);
dt = max(diff(t1))
t2 = t1(end)+dt:dt:p.Tmax;
t = [0 t1 t2];

t = [0 logspace(-4, log10(p.Tmax), 200)];
t(end) = p.Tmax;
% t = linspace(0, p.Tmax, 200);


tau = p.delta * t;
N = length(tau)

dt = 1e-7;
if (dt > min(diff(tau)))
    dt = min(diff(tau))
end

h = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);

hh = @(tau) -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + tau) / (p.beta - 1))))) * (p.beta - 1);
dh = @(t) imag(hh(t + 1i * 1e-6) / 1e-6);

% compute eigenvalues, eigenfunctions, adjoints
lambda = zeros(Nmodes, N);
lambda_a = lambda;
v = zeros(p.N, Nmodes, N);
adj = v;
for i = 1:N
    
    if asy
        
        for n = 0:Nmodes-1
           if (n == 0)
               v(:,n+1,i) = 1 / sqrt(h(i));
           else
               z = linspace(0, h(i), p.N);
               v(:,n+1,i) = sqrt(2 / h(i)) * cos(n * pi * z / h(i));
           end
        end
        
        adj(:,:,i) = v(:,:,i);
        
    else
        [ev, k, M] = comp_eigs(k0, p, h(i));
        
        [ef, ev] = eigs(M, Nmodes, 'SM');
        [tmp, P] = sort(diag(ev), 'descend');
        lambda(:,i) = tmp;
        v(:,:,i) = ef(:,P);
        
        Ma = comp_adj(k0, p, h(i));
        
        [ef, ev] = eigs(Ma, Nmodes, 'SM');
        [tmp, P] = sort(diag(ev), 'descend');
        lambda_a(:,i) = tmp;
        adj(:,:,i) = ef(:,P);
        
        
        % normalize eigenfunctions
        for j = 1:Nmodes
            v(:,j,i) = v(:,j,i) / sqrt(ip(v(:,j,i), v(:,j,i), h(i), p));
            adj(:,j,i) = adj(:,j,i) / ip(adj(:,j,i), v(:,j,i), h(i), p);
        end
    end
    
end

% go through and correct sign flips in the eigenfunctions
for i = 2:N
    for j = 1:Nmodes
        if (v(1,j,i) * v(1,j,i-1) < 0)
            v(:,j,i) = -v(:,j,i);
            adj(:,j,i) = -adj(:,j,i);
        end
        
    end
end

% asymptotic eigenvalues
if asy
    kh = k0 / sqrt(p.delta);
    C0 = -(1 - h) / p.beta ./ h;
    
    tmp = -kh^2 + ((h - 1 + p.beta) / p.beta) * kh^2 * p.Ma_prime / 80 - ...
        (1 - 2 * p.beta - 2 * p.beta * (1 - p.beta) * C0) ./ h;
    
    lambda(1,:) = p.delta * tmp;
    
    for k = 1:(Nmodes-1)
        tmp = -kh^2 - (1 - 2 * p.beta - 2 * p.beta * (1 - p.beta) * C0) .* 2 ./ h - ...
            (1 / 2 / pi^4 / k^4) * kh^2 * p.Ma_prime * (h - 1 + p.beta) / p.beta * (6 * (-1)^(k+1) + k^2 * pi^2 - 18);
        lambda(k+1,:) = -k^2 * pi^2 ./ h.^2 + p.delta * tmp;
    end
end



% interpolate eigenvectors for easy differentiation
if Nmodes > 1
    v_i = @(t_i) interp3(1:Nmodes, 1:p.N, tau, v, (1:Nmodes), (1:p.N), t_i);
else
    vv = reshape(v, p.N, N);
    v_i = @(t_i) interp2(tau, 1:p.N, vv, t_i, 1:p.N);
end


% plot the eigenvalues
if figs
    figure
    subplot(2,1,1);
    plot(tau, lambda);
    xlabel('\tau');
    ylabel('Eigenvalues');
    subplot(2,1,2);
    plot(tau, lambda_a);
    xlabel('tau');
    ylabel('Adjoint eigenvalues');
end


% plot some inner products
if figs
    figure
    for j = 1:Nmodes
        for i = 1:N
            tmp(i) = ip(v(:,j,i), v(:, j, i), h(i), p);
        end
        subplot(Nmodes, 1, j);
        plot(tau, tmp);
        xlabel('\tau');
        ylabel(sprintf('<v_%d, v_%d>',j,j));
    end
    
    figure
    for j = 1:Nmodes
        for k = 1:Nmodes
            for i = 1:N
                tmp(i) = ip(adj(:,j,i), v(:, k, i), h(i), p);
            end
            subplot(Nmodes, Nmodes, j + (k-1) * Nmodes);
            plot(tau, tmp);
            xlabel('\tau');
            ylabel(sprintf('<v_%d^*, v_%d>',j,k));
        end
    end
    
    figure
    for j = 1:Nmodes
        for i = 1:N
            tmp(i) = max(abs(adj(:,j,i)));
        end
        subplot(Nmodes, 1, j);
        plot(tau, tmp);
        xlabel('\tau');
        ylabel(sprintf('||v_%d^*||_{inf}',j));
    end
end

% compute numerical solution
% c0 = reshape(v(:,1,1), p.N, 1);
z = linspace(0, h(1), p.N)';
c0 = 16 * z.^2 .* (1 - z).^2;
num = stab(k0, p, c0);

% extract coeffs
num_i = deval(num, t);
c_num = zeros(Nmodes, N);
for i = 1:N
    for j = 1:Nmodes
        c_num(j,i) = ip(reshape(adj(:,j,i), p.N, 1), num_i(:,i), h(i), p);
    end
end

% compute the c_i
c_i = zeros(Nmodes, 1);
for i = 1:Nmodes
    c_i(i) = ip(reshape(adj(:,i,1), p.N, 1), c0, 1, p);
end
fprintf('c_i:\n');
disp(c_i);


% loop again to compute gamma
gamma = zeros(Nmodes, Nmodes, N);
gamma_asy = gamma;
dv = zeros(p.N, Nmodes);

% small-k gammas
d = dh(tau) ./ hh(tau);
for j = 0:Nmodes-1
    for k = 0:Nmodes-1
        
        if (j == 0 && k == 0)
            gamma_asy(j+1,k+1,:) = d / 2;
        elseif (j == k)
            gamma_asy(j+1,k+1,:) = d;
        elseif (k == 0)
            Ma = p.Ma_prime;
            kh = k0 / sqrt(p.delta);
            delta = p.delta;
            beta = p.beta;
            tmp = -delta * ((0.15e2 * h .^ 6 * kh ^ 2 * Ma * j ^ 2 * pi ^ 2 + kh ^ 2 * Ma * (j ^ 4 * pi ^ 4 + 0.120e3 + 0.20e2 * beta * j ^ 2 * pi ^ 2) * h .^ 5 + 0.80e2 * j ^ 4 * pi ^ 4 * beta) * ((-1) ^ (0.1e1 + j)) + (0.5e1 / 0.4e1 * kh ^ 2 * Ma * (j ^ 4 * pi ^ 4 + 0.24e2) * h .^ 6 + kh ^ 2 * Ma * (0.20e2 * j ^ 2 * pi ^ 2 + 0.120e3 * beta + j ^ 4 * pi ^ 4 * beta) * h .^ 5 - h .^ 3 * kh ^ 2 * Ma * j ^ 4 * pi ^ 4 / 0.4e1 + 0.20e2 * h * j ^ 4 * pi ^ 4 * beta + 0.80e2 * beta ^ 2 * j ^ 4 * pi ^ 4) * ((-1) ^ j) - 0.30e2 * kh ^ 2 * Ma * h .^ 6) .* sqrt(0.1e1 ./ h) * sqrt(0.2e1) / j ^ 6 / pi ^ 6 / beta / 0.20e2;
            gamma_asy(j+1,k+1,:) = tmp .* dh(tau);
        elseif (j == 0)
            gamma_asy(j+1,k+1,:) = sqrt(2) * (-1)^(k) * d;
        else
            gamma_asy(j+1,k+1,:) = 2 * (-1)^(j+k+1) * k^2 / (j^2 - k^2) * d;
        end
    end    
end

if (asy)
    gamma = gamma_asy;
else
    for i = 1:N
        
        
        h_2 = interp1(tau, h, tau(i) + dt);
        v_2 = v_i(tau(i) + dt);
        h_1 = interp1(tau, h, tau(i) - dt);
        v_1 = v_i(tau(i) - dt);
        
        for k = 1:Nmodes
            
            if (i == 1)
                tmp_2 = interp1(linspace(0, h_2, p.N)', v_2(:,k), linspace(0, h(i), p.N)', 'pchip','extrap');
                dv(:,k) = (tmp_2 - v(:,k,i)) / dt;
            elseif (i == N)
                tmp_1 = interp1(linspace(0, h_1, p.N)', v_1(:,k), linspace(0, h(i), p.N)', 'pchip','extrap');
                dv(:,k) = (v(:,k,i) - tmp_1) / dt;
            else
                tmp_2 = interp1(linspace(0, h_2, p.N)', v_2(:,k), linspace(0, h(i), p.N)', 'pchip', 'extrap');
                tmp_1 = interp1(linspace(0, h_1, p.N)', v_1(:,k), linspace(0, h(i), p.N)', 'pchip', 'extrap');
                dv(:,k) = (tmp_2 - tmp_1) / 2 / dt;
            end
            
            if (max(isnan(dv(:,k))) == 1)
                fprintf('NaN detected\n');
                fprintf('i = %d, h_1 = %.4e, h_2 = %.4e, v_1 = %d, v_2 = %d\n', i, h_1, h_2, max(max(isnan(v_1))), max(max(isnan(v_2))));
                return
            end
        end
        

        
        for j = 1:Nmodes
            for k = 1:Nmodes
                gamma(j,k,i) = -ip(reshape(adj(:,j,i), p.N, 1), dv(:,k), h(i), p);
            end
        end
        
        
    end
end

% compute the c by direct numerical simulation
c_ns = c_solve(c_i, lambda, gamma, tau, p);


% plot gamma
if figs
    figure
    for j = 1:Nmodes
        for k = 1:Nmodes
            subplot(Nmodes, Nmodes, k + (j - 1) * Nmodes);
            tmp = [reshape(gamma(j,k,:), N, 1) reshape(gamma_asy(j,k,:), N, 1)];
            plot(tau, tmp);
            xlabel('\tau');
            ylabel(sprintf('\\gamma_{%d%d}',j,k));
        end
    end
end

% compute beta
tmp = zeros(N, Nmodes);
for i = 1:Nmodes
    tmp(:,i) = reshape(gamma(i,i,:), N, 1);
end
beta = cumtrapz(tau', lambda' / p.delta + tmp)';

% plot beta
if figs
    figure
    plot(tau, beta);
    xlabel('\tau');
    ylabel('\beta');
end

% compute leading order coeffs
c_lo = zeros(Nmodes, N);
for i = 1:Nmodes 
    c_lo(i,:) = c_i(i) * exp(beta(i,:));
end

fprintf('c_lo(t = 0):\n');
disp(c_lo(:,1));
fprintf('c_num(t = 0):\n');
disp(c_num(:,1));

% do the corrections
c_corr = zeros(Nmodes, N);
for j = 1:Nmodes
    for l = 1:Nmodes
        if (l ~= j)
            int = cumtrapz(tau, reshape(gamma(j, l,:) .* gamma(l, j,:), 1, N) ./ (lambda(j,:) - lambda(l,:)));
            
            tmp = -c_i(l) * reshape(gamma(j, l, :), 1, N) ./ (lambda(j, :) - lambda(l, :)) .* exp(beta(l,:)) + ...
                (c_i(j) * int + gamma(j, l, 1) * c_i(l)  / (lambda(j,1) - lambda(l,1))) .* exp(beta(j,:));
                        
            c_corr(j,:) = c_corr(j,:) + tmp;
        end
    end
end

% t_i = logspace(log10(c_ns.x(2)), log10(c_ns.x(end)), 30);
t_i = linspace(c_ns.x(2), c_ns.x(end), 30);

% plot the amplitudes
A_num = max(abs(num.y));
A_lo = zeros(1,N);
A_corr = A_lo;
for i = 1:N
    A_lo(i) = max(abs( reshape(v(:,:,i), p.N, Nmodes) * c_lo(:,i)  ));
    A_corr(i) = max(abs( reshape(v(:,:,i), p.N, Nmodes) * (c_lo(:,i) + p.delta * c_corr(:,i)) ));
    A_ns(i) = max(abs(reshape(v(:,:,i), p.N, Nmodes) * deval(c_ns, t(i))));
end
if figs
    figure
    plot(num.x * p.delta, A_num / norm(c0, 'inf'), 'k', tau, A_lo / norm(c0, 'inf'), 'b-.', tau, A_corr / norm(c0, 'inf'), 'r--');
    hold on;
    plot(t_i * p.delta, interp1(tau, A_ns, t_i * p.delta) / norm(c0, 'inf'), 'mo');
    xlabel('\tau');
    ylabel('Amplitudes');
    l = legend('numerical','leading order','leading plus correction', 'num soln ode system', 'location', 'best');
    figure
    loglog(num.x * p.delta, A_num / norm(c0, 'inf'), 'k', tau, A_lo / norm(c0, 'inf'), 'b-.', tau, A_corr / norm(c0, 'inf'), 'r--');
    xlabel('\tau');
    ylabel('Amplitudes');
    l = legend('numerical','leading order','leading plus correction', 'location', 'best');
end

% rescale amps
A_num = A_num / norm(c0, 'inf');
A_lo = A_lo / norm(c0, 'inf');
A_corr = A_corr / norm(c0, 'inf');
A_ns = A_ns / norm(c0, 'inf');


% plot the coefficients
if figs
    figure
    for j = 1:Nmodes
        subplot(ceil(Nmodes/2),2,j);
        semilogx(tau, c_num(j,:), 'k', tau, c_lo(j,:), 'b-.', tau, c_lo(j,:) + p.delta * c_corr(j,:), 'r--');
        hold on;
        semilogx(t_i * p.delta, deval(c_ns, t_i, j), 'mo');
        xlabel('\tau');
        ylabel(strcat('c_',num2str(j)));
        tmp = get(gca, 'xlim');
        xlim([tmp(1), p.Tmax * p.delta]);
    end
end

% log-plot the coefficients
if figs
    figure
    for j = 1:Nmodes
        subplot(ceil(Nmodes/2),2,j);
        loglog(tau, abs(c_num(j,:)), 'k', tau, abs(c_lo(j,:)), 'b-.', tau, abs(c_lo(j,:) + p.delta * c_corr(j,:)), 'r--');
        hold on;
        loglog(t_i * p.delta, abs(deval(c_ns, t_i, j)), 'mo');
        xlabel('\tau');
        ylabel(strcat('c_',num2str(j)));
        tmp = get(gca, 'xlim');
        xlim([tmp(1), p.Tmax * p.delta]);
    end
end


% compare coeffs from ode solve and from asymptotics
if figs
    figure
    for j = 1:Nmodes
        subplot(ceil(Nmodes/2),2,j);
        loglog(tau, abs(c_lo(j,:) + p.delta * c_corr(j,:) - deval(c_ns, t, j)));
        xlabel('\tau');
        ylabel(sprintf('|c_%d - c_{%d,ode}|',j,j));
        tmp = get(gca, 'xlim');
        xlim([tmp(1), p.Tmax * p.delta]);
    end
    
    % compare coeffs from ode solve and from full numerics
    figure
    for j = 1:Nmodes
        subplot(ceil(Nmodes/2),2,j);
        loglog(tau, abs(c_num(j,:) - deval(c_ns, t, j)));
        xlabel('\tau');
        ylabel(sprintf('|c_{%d, num} - c_{%d,ode}|',j,j));
        tmp = get(gca, 'xlim');
        xlim([tmp(1), p.Tmax * p.delta]);
    end
end